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This paper describes a computer analysis of the propagating modes 
of a rectangular dielectric waveguide. The analysis is based on an expansion 
of the electromagnetic field in terms of a series of circular harmonics, that is, 
Bessel and modified Bessel functions multiplied by trigonometric functions. 
The electric and magnetic fields inside the waveguide core are matched to 
those outside the core at appropriate points on the boundary to yield equa- 
tions which are then solved on a computer for the propagation constants and 
field configurations of the various modes. 

The paper presents the results of the computations in the form of curves of 
the propagation constants and' as computer generated mode patterns. The 
propagation curves are presented in a form which makes them refractive- 
index independent as long as the difference of the index of the core and the 
surrounding medium is small, the case which applies to integrated optics. 
In addition to those for small index difference, it also gives results for 
larger index differences such as might be encountered for microwave appli- 
cations. 

I. INTRODUCTION 

It is anticipated that dielectric waveguides will be used as the 
fundamental building blocks of integrated optical circuits. These wave- 
guides can serve not only as a transmission medium to confine and 
direct optical signals, but also as the basis for circuits such as filters 
and directional couplers. 1 Thus, it is important to have a thorough 
knowledge of the properties of their modes. 

Circular dielectric waveguides have received considerable attention 
because circular geometry is commonly used in fiber optics. 2-8 In many 
integrated optics applications it is expected that waveguides will con- 
sist of a rectangular, or near rectangular, dielectric core embedded in 
a dielectric medium of slightly lower refractive index. The modes 
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for this geometry are more difficult to analyze than those of the me- 
tallic rectangular waveguide because of the nature of the boundary. 

Marcatili, using approximations based on the assumption that most 
of the power flow is confined to the waveguide core, has derived in 
closed form the properties of a rectangular dielectric waveguide. 6 In 
his solution, fields with sinusoidal variation in the core are matched 
to exponentially decaying fields in the external medium. In each 
region only a single mode is used. The results of this method are 
obtained in a relatively simple form for numerical evaluation. 

The properties of the principal mode of the rectangular dielectric 
waveguide have been studied by Schlosser and Unger using a high- 
speed digital computer. 7 In their approach the transverse plane was 
divided into regions, as shown in Fig. 1, and rectangular coordinate 
solutions assumed in each of the regions. The longitudinal propagation 
constant was then adjusted so that a field match could be achieved 
at discrete points along the boundary. This method gives results 
which, theoretically, are valid over a wider range than Marcatili's, 
but with a significant increase in computational difficulty. One short- 
coming of the method is that for a given mode, as the wavelength 
increases the field extent increases, so, in the limit it becomes increas- 
ingly difficult to match the fields along the boundaries between regions 
[1] and [2] and between regions [2] and [3]. 

A variational approach has been undertaken by Shaw and others. 8 
They assume a test solution with two or three variable parameters 
in the core. From this test solution, the fields outside the core are 
then derived and the parameters are varied to achieve a consistent 
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Fig. 1 — Matching boundaries for rectangular mode analysis. 
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solution. This approach, like that of Schlosser, requires involved com- 
putations. Also, it has the disadvantage that the test function must 
be assumed in advance. In addition, some of his preliminary results 
do not show the proper behavior for the limiting cases (waveguide 
dimensions which are very large or very small compared with the 
wavelength) . 

In the present analysis the radial variation of the longitudinal 
electric and magnetic fields of the modes are represented by a sum 
of Bessel functions inside the waveguide core and by a sum of modi- 
fied Bessel functions outside the waveguide core. Solutions are found 
by matching the fields along the perimeter of the core. Thus, the 
matching boundary is not a function of the waveguide parameters, 
so the computational complexity does not increase with wavelength. 

Section II discusses the underlying theory of the circular-harmonic 
analysis of rectangular dielectric waveguides. This is followed by a 
description of computational techniques and special graphical methods 
of presentation used. Section III is divided into three parts, the first 
describing the accuracy of the computations, the second describing 
field patterns, and the third presenting propagation curves. 

II. DERIVATION OF EQUATIONS 

The waveguide considered here consists of a rectangular core of 
dielectric constant, ei, surrounded by an infinite medium of dielectric 
constant, e . Both media are assumed to be isotropic, and have the 
permeability of free space, /*<>. Figure 2 shows the coordinate systems 
(rectangular and cylindrical) and rod dimension used in this paper. 
The direction of propagation is iv the + z direction (towards the 
observer) . 

In cylindrical coordinates the field solutions of Maxwell's equations 
take the form of Bessel functions and modified Bessel functions mul- 
tiplied by trigometric functions, and their derivatives. In order for 
propagation to take place in the z direction, the field solutions must 
be Bessel functions in the core and modified Bessel functions outside. 
Since Bessel functions of the second kind have a pole at the origin 
and modified Bessel functions of the first kind a pole at infinity, the 
radial variation of the fields is assumed to be a sum of Bessel func- 
tions of the first kind and their derivatives inside the core and a sum 
of modified Bessel functions and their derivatives outside the core. 

In cylindrical coordinates, the z components of the electric and 
magnetic fields are given by 
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E tl = ^ o-nJnQir) sin (nd + #„) exp [i(k t z — at)] 



and 



H tl = ^ b„J„(hr) sin (nd + &■) exp [t(&,z — w/)] 

n-0 

inside the core, and by 

co 

E,o = ^ c n K n (pr) sin (n0 + ^ n ) exp [t'(^«z — «0] 

n-0 

and 

00 

H t0 = ^2 d n K n (pr) sin (nd + ^„) exp [t'(/c,3 — wt)] 



(la) 



(lb) 



(lc) 



(Id) 



outside the core, where u> is the radian frequency and k t the longitudinal 
propagation constant. The transverse propagation constants are given 
by 



h=(kl- k 2 ) h 



(2a) 



and 



V - <K - fcoV (2b) 

where /^ = &>Gmi)* an d ^o = coGuoCo) 1 - The terms «/„ and K n are the nth 
order Bessel functions and modified Bessel functions, respectively, and 
^ n and <p„ are arbitrary phase angles. 

The transverse components of the fields are given by 9 



_ %K [§E, Azg«\ dH.1 
*' ~ k 2 - L dr + \k,r) dd § ] 



(3a) 
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Fig. 2 — Dimensions and coordinate system. 
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(3b) 



1 dE t U^\ BH z 

r dd \ k, ) dr 



r - kt 



( k 2 \dE z dH,~] ,,, 

-l^J^ + irj (3c) 

( k 2 \ dE z 1 BRA ,,« 



where k can be either fc x or fc . 

Finally, the component of the electric field tangent to the rectang- 
ular core is given by 

E t = ±(E T sm 6 + E 6 cos 6) ~° e < ° < ° e (4a) 

7T - 6 C < <7T + e r 

or 



E, = ±(-E r cos 6 + E g sin 6) , 6 C < d < tt - 

t + e e < e < - e e 



(4b) 



where 6 C is the angle which a radial line to the corner in the first 
quadrant makes with the x axis. Similar expressions exist for the 
tangential magnetic field. 

2.1 Effects of Symmetry 

Since the waveguide is symmetrical about the x axis the fields 
must be either symmetric or antisymmetric about this axis. This is 
true because the structure is invarient under 180° rotations and there- 
fore the field patterns must be invarient under a 180° rotation, except 
for sign. From this and the fact that d/dd appears in each of equations 
(3), it is evident that two types of modes must exist, the first type 
with <p„ = and \p n = tt/2 and the second type with <p n = ir/2 and 

Similarly, the field functions must also be symmetric or anti- 
symmetric about the y axis. Suppose, for example, E~o exhibits a sinu- 
soidal angular dependence about 6 = {E :n is odd about the x axis) . 
Then, letting a = 6 — 7r/2, equation (lc) can be put in the form 

00 

E, = ^ c n K n (pr)(s'wna cosnr/2 + cos na sin nr/2) . (5) 

n-0 

For E g0 to be purely symmetric about a = (the y axis), all n must 
be odd; for E~ to be antisymmetric about a = all n must be even. 
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Since similar results apply for cosinusoidal variation of E^ about 
= 0, and all other field functions as well, any given mode must 
consist of either even harmonics or odd harmonics. 

From the preceeding analysis it is evident that if the matching 
points are selected symmetrically about both the x and y axes, then, 
except possibly for sign, every point will have an equivalent point 
in each quadrant. Therefore, the field matching need only be per- 
formed in one quadrant. Thus, the use of the symmetry of the struc- 
ture not only reduces the number of constants required to calculate 
the properties of a given mode by a factor of four, it also decreases 
the number of points to achieve a given degree of accuracy by the 
same factor. 

2.2 Selection of Matching Points 

As mentioned in Section 2.1, the matching point locations should be 
symmetrical about the x and y axes. For the odd harmonic cases, the 
points used to compute the results to be presented in Section III were 
6 m = (m — l/2)ir/2N; m — 1, ■ • • , N, where N was the number of 
space harmonics. 

The choice of points for the even harmonic cases was more complicated 
since simultaneous existence of an n = harmonic for both the TE and 
TM circular modes is inconsistent with the waveguide symmetries. 
Thus, if the maximum n for both the TE and TM solutions are equal, 
the total number of coefficients to be found will be 4N — 2 rather than 
4=N as in the previous case. 

The method of choosing points for the even harmonic modes used for 
the computation of the results of Section III was to pick the points 
for the field components with even symmetry about = to be m = 
(m — l/2)ir/2N; m = 1, 2, • • • , N, and for the field components with 
odd symmetry about = to be m = (m — N — l/2)ir/2(N — 1); 
m = N + 1, N + 2, ■ ■ ■ , (2N — 1) for cases with unity aspect ratio, 
(a/b = 1). For aspect ratios other than unity, all points were chosen 
according to the first formula, except that the first and last points for 
the odd z component were omitted. 

2.3 Formulation of Matrix Elements 

The coefficients of equation (1) were found by matching the tan- 
gential electric and magnetic fields along the boundary of the wave- 
guide core. Since each type of field consists of both longitudinal and 
transverse components, four types of matching equations exist. 

To facilitate computer analysis the matching equations were put in 
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matrix form. The matching equations in matrix form for the longi- 
tudinal field components are 

E LA A = E LC C (6a) 
for the electric field and 

H LD B = H LP D (6b) 

for the magnetic field. For the transverse fields the matrix matching 
equations arc given by 

E TA A + E T "B = E TC C + E Tn D (6c) 
for the electric field and 

H TA A + H TB B = H TC C + H TD D (6d) 

for the magnetic field. The .4, B, C, and D matrices are N element 
column matrices of the a n , b„, c„, and d n mode coefficients, respectively. 
The elements of the m X n matrices E LA , E LC , H LB , H lD , E TA , E Ta , 
E TC , E TD , H TA , H TB , H TC , and H TD are given by 

ett = JS, (7a) 

ett = KS, (7b) 

ht" = JC, (7c) 

*£? = KC, (7d) 

elt = -UJ'SR + JCT), (7e) 

tS - fc n Zo(J,Si? + J'CT), (7f) 

rl n c = fc,(K'£i2 + KCr), (7g) 

e™ = -k Z (KSR + K'CT 1 ), (7h) 

fc" = eMK'R - J'ST)/Z n , (7i) 

h™ = -KiJ'CR - J ST), (7j) 

hZ = -k (KCR - K'ST)/Z , (7k) 

hl D = fc.(K'CB - KST), (71) 



.'here 



£ = (Mo/co) 4 , 
e r = e,/e» , 
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5 = sin (nd m + <p) ( tp = 
| or 

C = cos (nd m + <p) tp = tt/2 

./ = J„(/ir m ), A' = A-„(pr m ), 



and 





./' = 


= J'n(hr m ), 


K' = K' n (pr, n ), 




T n,J n (hr m ) 


nK„(pr„) 

&. — 2 
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T , J!,(hr m ) 
J h ' 


K/ Ki(pr m ) 
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R = sin m 






R = —cos 6 m 


T = cos 0„ 




\ d < e e , 


T = sin d m 


r. = (a/2) 


cos 6 m 




r m = (6/2) sin 



> 



For 6 = d e , the boundary at the corner was assumed to be perpendicular 
to the radial line connecting it to the origin, so for this case R = cos (0 m 
+ tt/4), T = cos (d m - t/4), and r m = (a 2 + b')*/4. 

2.4 Mode Designation 

Unlike metallic waveguides, the field patterns of dielectric wave- 
guides are sensitive to refractive index difference, wavelength, and 
aspect ratio. This complicates the problem of finding a reasonably 
descriptive mode designation scheme. 

For rectangular metallic waveguides, the accepted approach is to 
designate the modes as TE (or H) and TM (or E), and to specify 
the number of field maxima in the x and y directions with a double 
subscript. When there is no variation the subscript is used. 

Since the rectangular dielectric waveguide modes are neither pure 
TE nor pure TM, a different scheme must be used. The scheme adopted 
is based on the fact that in the limit, for large aspect ratio, short wave- 
length, and small refractive index difference, the transverse electric 
field is primarily parallel to one of the transverse axes. Modes are 
designated as E y mn if in the limit their electric field is parallel to the y 
axis and as E z mn if in the limit their electric field is parallel to the x axis. 
The m and n subscript are used to designate the number of maxima 
in the x and y directions, respectively. 1. 



t This scheme agrees with that used by Marcatili in Ref. 6. 
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2.5 Electric and Magnetic Field Function Differences 

For a hollow metallic waveguide where pure TE and TM modes 
can exist, it is evident from equation (3) that E r and He have similar 
transverse variations as do Eq and H r , so that the impedance is in- 
dependent of position. Furthermore, the transverse electric and mag- 
netic fields are perpendicular and the power flow, Re {E X H*}, does 
not change sign anywhere across the waveguide. 

By examination of equation (3), it is clear that for the mixed modes 
of the dielectric waveguide, the field functions are not similar and the 
impedance is a function of position. In order for the transverse fields 
E t and H t to be perpendicular, 

E,H, = E T Hr + E e H s = 0. (8) 

Now, from equation (3) 

k\ - k 2 (dH, dE. 1 dH, BE\ 
E%Et ~ k\ \dr dr + r 2 dd dd ) W 

Thus, E t and H t are not necessarily perpendicular. Finally, since the 
transverse variations of E t and H t are not the same, the electric field 
and magnetic field can change sign at different points, which results in 
negative power flow.* 

Three special cases exist where the electric and magnetic fields, and 
the impedance, have the same positional dependence, and where the 
power flow does not change sign across the waveguide: 

(i) in one of the regions if the propagation constant is approximately 
equal to the bulk propagation constant of that region, that is, if k tt k x 
or k tt k , 

(ii) everywhere in the limit for small refractive index difference, 
since case i will then hold in both regions, and 

(Hi) everywhere for circular symmetry of both the structure and the 
modes. 

2.6 Normalization 

The arguments of the Bessel and modified Bessel functions are given 
by hr = (fc 2 — fc 2 )V and pr = (k 2 t — hfflr, respectively. The first argu- 
ment can be put in the form 

hr = [hi -k 2 - p 2 ]V (10) 



t This unusual property has also been observed for helices. 10 Presumably, if 
loss were included there would be a radial component of power to feed the re- 
verse flow, and the lossless case can be thought of as the limit of the lossy case. 
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Letting 

(P 2 = (fe ffi] ~ * , (11) 

and 

(ft = >* (n 2 ~ 1) J , (12) 

where 

rc r = (fe./fco)* (13) 

is the index of refraction of the core relative to the outer medium, gives 

?)/• = (P(R (14) 

and 

hr = (R(l - (P 2 )*. (15) 

The curves of the propagation constant given in Section III are 
drawn in terms of (P 2 and (B, where 

CB = ^ (n 2 . - 1)* (16) 

Ao 

and X = 27r/fc . Since (ft is proportional to 1/(71?. — 1) and (P and (B 
are proportional to (n 2 . — 1)*, the use of (P 2 and (B as plotting variables 
eliminates the explicit dependence of the Bessel and modified Bessel 
function arguments on the refractive indices of the media. 

Examination of the matching equations, equations (6), reveals that 
e r appears in the H TA term. However, since e r appears as a multiplicative 
factor in H TA , for sufficiently small values the normalized propagation 
constant, (P 2 , is independent of e r . 

The normalized propagation constant, (P 2 , has two additional prop- 
erties which make its use convenient. First, its range of variation is on 
the interval (0, 1). Second, for n T tt 1, 

(P 2 « k ' /ko - 2 , (17) 

An T 

where An r = n r — 1; so for small n r , (P 2 is proportional to k t — ko . The 
latter property is the reason that <P 2 rather than (P was used as a plotting 
variable. 

2.7 Method of Computation 

2.7.1 Propagation Constant 

Equation (6) yields 4N simultaneous homogeneous linear equations 
for the a n) b„, c„, and d n for the odd modes and 4iV"-2 equations for 
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the even modes, using the matching points previously described. The 
equations can be combined to form a single matrix equation 



[Q][T\ = 0, 



(18) 



where 



Q = 



and the column matrix 



/<:' 



o 



-E l 



H LB 

E TA E T " -E TC 

H TA H TB __ H TC 




-H LD 

-E TD 

_ H TD 



[T] = 



All of the quantities in the matrices [Q] and [T] are themselves ma- 
trices as defined by equations (1) , (6) , and (7) . 

In order for a nontrivial solution to equation (18) to exist 



Det [Q] = 0. 



(19) 



The normalized propagation constant, (P 2 , was found by substituting 
test values into equation (19). First, values of (P 2 evenly distributed in 
the interval (0, 1) were substituted to crudely locate the roots. Then, 
Newton's method was used to find the roots to the desired accuracy. 11 
Generally, one Newton approximation was used to find (P 2 for the prop- 
agation curves and about ten Newton's approximations when (P 2 was 
to be used to calculate field plots. 

Both the simple method of triangulation 12 and the more complicated 
Gauss pivotal condensation method 13 were used to evaluate the deter- 
minant, the former for almost all cases and the latter for a few cases 
when roundoff error was apparent because the value of the determinant 
was not a smooth function of (P 2 . In all cases double precision arithmetic 
was used. For five space harmonics, about 0.1 second of IBM 360/65 
computing time was required for each value of (P 2 to evaluate the deter- 
minant using the triangulation method. 

Due to the wide dynamic range of the coefficients, steps had to be 
taken to prevent underflow and overflow of the computer and to re- 
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duce the effects of roundoff. Multiplying a row or column of the ma- 
trix by a finite constant is equivalent to multiplying the determinant 
by that constant. Thus, any row or column of the determinant can be 
multiplied by a positive function without shifting its zeroes. 

A detailed theory giving the "best functions" can be derived. How- 
ever, since a "brute force" method was used, the more sophisticated 
method, which was not used because it would have required a substan- 
tial increase in the complexity of the program logic, is not discussed. 
It was found that multiplying the Bessel function terms by h 2 d/\J n (hb) \ 
and the modified Bessel function terms by p 2 d/k n (pb), where d is 
the average of the waveguide dimensions, kept the variation of the 
terms "under control." A further simplification was made by setting 
Z to unity, which does not shift the zeroes of the determinant be- 
cause if the H t rows are multiplied by Z , then if Z appears in a 
column, it will appear in a similar manner in every element of the 
column. 

2.7.2 Mode Configurations 

The electric and magnetic fields were calculated for representative 
cases from equation (3). To find the a.,,, b„, c n , and d n coefficients, 
k z was first found from equation (19). Its value was then substituted 
into equation (18). By setting one of the elements of the T column 
matrix to unity, all of the other elements were then found by standard 
matrix techniques. 13 

Several approaches were used to obtain information that could 
be used to derive the field patterns. These included computation of 
the field components along radial cuts of the waveguide cross section, 
computer generated isoclines giving the direction of the electric field, 
and computer generated mode pictures. 

The isoclines and pictures were drawn using a simulated Stromberg 
Carlson SC-4020 cathode ray tube plotter, which is capable of gen- 
erating points and lines on a 1024 X 1024 grid.* A single quadrant 
was used for the isoclines and intensity picture since the results for 
all quadrants are identical except for orientation. In general, the di- 
mensions were scaled so that the long dimension of the rectangular 
waveguide core extended over 80 percent of the displayed width. All 
figures were plotted at the points (20m, 20n) , where m and n take on 
all integer values from to 49. 

Isocline drawings were made by drawing a line at each of the co- 
ordinate points parallel to the electric field at that point (all lines 



t An SC-4060 plotter was used to simulate the SC-4020 plotter to take advan- 
tage of previously existing programs. 
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had the same length). The isocline drawings were used as working 
tools to derive the field line drawings in Section III. 

In order to draw pictures of mode patterns, the power density was 
calculated at each of the points to be plotted. The square root of the 
power density was then normalized to the square root of the peak 
power density and quantized into 21 levels. About each point in the 
picture, a portion of the figure shown in Fig. 3 was then plotted, 
starting at 1 and going to the point corresponding to the appropriate 
quantized level (except at the points where the quantized power 
was zero where no plotting was done) . Since the size of the cathode 
ray tube spot is approximately equal to the line spacing in the figure, 
the plotted figures are filled in. Therefore, the light passed by these 
figures is approximately equal to the power density to be represented. 
For small index difference, since the power density is proportional to 
the square of the transverse electric field, the dynamic range of the 
pictures (in terms of the electric field) is 400. 

Starting with the single quadrant pictures, complete pictures were 
generated by making quadruple exposures of the microfilm. In general, 
about 30 to 60 seconds of IBM 360/65 computing time were required 
for each picture. 

III. RESULTS OP COMPUTATION 

This section gives the computed results. Section 3.1 discusses ac- 
curacy. This is followed by a discussion of field plots and mode 
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Fig. 3 — Intensity picture figure. 
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Table I — Sample Accuracy Results 




Number of 


(p2 




a/b = l 


a/6 =2 


a/b=3 


a/6 =4 


3 
4 
5 
6 
7 
8 
9 


0.714 
0.713 
0.715 
0.714 
0.715 
0.715 
0.715 


0.811 
0.811 
0.808 
0.808 
0.808 
0.807 
0.807 


0.820 
0.820 
0.819 
0.822 
0.820 
0.820 
0.823 


0.828 
0.819 
0.813 
0.820 
0.813 
0.814 
0.815 


Variation 


0-2% 


0.4% 


0.4% 


1.5% 



pictures in Section 3.2. Finally, curves of the propagation constant for 
a variety of conditions are presented in Section 3.3. 

3.1 Accuracy 

Numerous test runs were made in order to obtain an estimate of 
the accuracy of the computed results. The results of several of these 
runs are given in Table I for the first mode with (B = 2. The numbers 
at the bottom of the table represent the total variation for a given aspect 
ratio taken as a percentage of the full range possible (one). 

For small aspect ratios, it is clear that the convergence is very rapid. 
However, for larger aspect ratios the convergence is not as good. For 
example, the variation for an aspect ratio of four is 1.5 percent (taken 
as a percentage of the full range of variation). For this case, from the 
table and from the limit for infinite aspect ratio 14 which is an upper 
bound for (P 2 , it appears the error is about 3 percent. This error is 
achieved with a relatively small number of harmonics and can only 
be improved by using a prohibitively large number of harmonics on a 
computer which carries more significant digits than the one which was 
available for this study. However, since solutions exist for an infinite 
aspect ratio, the decrease in accuracy for the large aspect ratio of the 
circular-harmonic method is not a serious problem. 

Computations similar to those for Table I were performed to ob- 
tain an estimate of the upper bound of the accuracy of the cases pre- 
sented in Section 3.3. From these calculations, it is believed that all 
of the data to be presented in the following sections is accurate to 
1 percent, except for the results of calculations using even harmonics 
for aspect ratios other than unity which are believed to be accurate 
to better than 2 percent. In general, accuracy decreases as the mode 
order increases, although not monotonically. 
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The results of the circular-harmonic analysis and of Marcatili's 
analysis agree. In the regions where his method and the circular- 
harmonic method are both theoretically valid, the agreement is well 
within the tolerances given above. To avoid duplication, the reader is 
directed to his curves for a comparison. 

The effect of the number of harmonics used in the field patterns is of 
some interest. This question has not been explored in great detail; 
however, a few comparisons of intensity pictures for different numbers 
of circular harmonics were made. In general, it was found that five 
harmonics were sufficient to give a good representation of the modes 
that this paper presents. An example of this is given in Fig. 4, comparing 
the E\ x mode intensity patterns for five and nine harmonics. For the 
results which follow, five circular harmonics were used. 

3.2 Mode Configurations 

Figure 5 shows intensity pictures for the first six modes for unity 
aspect ratio, (B = 3, and an index difference of 0.01. Figure 6 gives 
similar data for an aspect ratio of two and (B = 2. For both, the plots 
are arranged in ascending order of cutoff frequency. All of the pictures 
are for E" mn modes. These pictures are virtually indistinguishable from 
the corresponding E" mn modes so both sets are not presented. In general, 
for small index differences the E v mn and E x mn can be considered to be 
near duals, that is, to have identical field patterns except that the 
electric and magnetic fields are interchanged. 

The field distribution patterns for the modes of Figs. 5 and 6 are 
more complicated than those for the rectangular metallic waveguide 




Fig. 4 — Intensity for the E'{ v mode for a/b = 2, flj = 2, and An r = .01 : (a) for 
five harmonics and (b) for nine harmonics. 
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Fig. 7 — Field configuration of the E u lx mode. 

since they extend beyond the waveguide boundary and, in general, their 
shape is dependent on waveguide parameters other than shape. The E x Xi 
and E\ x modes have the simplest field patterns. Figure 7 shows the elec- 
tric and magnetic field orientations for the E"n mode. In this figure and 
the following ones, there are heavy lines in the regions of high field inten- 
sity and light lines in regions of low field intensity. Only E v mn modes are 
shown since the E x mn modes can be obtained by interchanging the electric 
and magnetic field vectors. 

Figure 8 shows the field lines for the E v 2l and E\ 2 modes for a large 
aspect ratio. (For a/b — » «> the fields have the appearance of rectangular 
metallic waveguide modes.) However, as the aspect ratio approaches 
unity, the E\ 2 and E 21 modes and the E 21 aiQ d EU modes couple and 
shift to the patterns shown in Fig. 9. Most of the change takes place 
with the aspect ratio close to unity. 

Figures 10, 11, and 12 show the field configurations for the E 22 mode, 
the E\ x mode, and the E" a mode, respectively. The field patterns of 
these modes do not change drastically with the aspect ratios. 

Figure 13a shows an intensity picture of the E\ 2 mode and Figure 




+-:£££-+ 
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(a) (b) 

Fig. 8 — Field configurations for the (a) E^ and (b) E\ 2 modes far from cutoff. 
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Fig. 9 — Field configurations for the square (a) E% t and (b) E\ 2 modes. 

13b its field pattern for unity aspect ratio. The field pattern inside the 
core is similar to a sum of the TE 23 and TE 32 of metallic waveguide, 
shown in Fig. 13c and d, respectively. Figure 13a demonstrates that the 
circular-harmonic analysis can generate complex field patterns with 
a relatively small number of harmonics. 

Figures 14 and 15 show the variation of the intensity distribution with 
(P 2 for the E v xl and E"i modes, respectively. As one would expect, for 
small values of (P 2 the radial extent of both modes increases very rapidly 
as (P 2 decreases. It is of significance, however, that most of the energy is 
contained within the waveguide core, even for relatively small values 
of (P 2 and An. Thus, Marcatili's assumption that very little energy 
propagates in the region of the corners is valid over a wide range. 

3.3 Propagation Curves 

In all cases of computed propagation curves, the normalized wave- 
guide height (B, as given in equation (11), is plotted on the horizontal 




Fig. 10 — Field configuration of the E" 2i mode. 
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Fig. 11 — Field configuration of the £g, mode. 

axis and the normalized propagation constant, (P 2 , given in equation (16), 
along the vertical axis. 

Figure 16 shows the case of vanishing index difference for an aspect 
ratio of one. The first 16 modes are shown. For this case the following 
six degenerate groups ,exist 

■C'll , ■C'll 

JJTV Tpx 77TW Tpx 

•^12 » -"12 . ^21 » ^21 
-^31 » E* 3 

Ell , E\ 3 

EX T?V 

22 . -"22 

-"32 ) -"23 i -"23 i -"23 • 

In addition, the E\ x and the E* zy modes are almost degenerate except 




Fig. 12 — Field configuration of the E\ 3 mode. 



WAVEGUIDE ANALYSIS 



2153 







Fig. 13 — The E% 2 mode for unity aspect ratio: (a) intensity, (b) field configura- 
tion, (c) TE n , and (d) TE U . 

near cutoff. The splitting of these modes can be accounted for by the 
differences of the field patterns shown in Fig. 11 and 12. Since the E\ x 
mode reversals occur along the direction of the electric field lines, the 
electric field for this mode must have a larger longitudinal field com- 
ponent than for the E\ x mode. 

All degeneracies, except the E v nn — E' mn , are broken by a change 
in the aspect ratio as demonstrated in Fig. 17, which is drawn for 
the first 12 modes of a waveguide of aspect ratio 2. One interesting 
feature of this curve is the mode crossing of the E v 31 and E\ 2 modes. 
Crossings of this type, which cannot occur in metallic waveguides, are 
possible because the field functions are frequency dependent. Qualita- 
tively, it can be explained by noting that field reversals must take place 
in the core, therefore constraining the central lobe of the E v al more than 
any of the E\ 2 mode lobes as cutoff is approached. Far from cutoff, 
however, all fields are well constrained and the E^ t mode has a larger 
propagation constant than the E\ 2 mode, as it does for the similar 
metallic waveguide mode with an aspect ratio of 2. 

The effect of finite index difference on the modes can be observed by 
comparing Fig. 16, which is computed for unity aspect ratio and a 
vanishing index difference, with Fig. 18, which is computed for unity 
aspect ratio and a 0.5 index difference. The curves for modes whose 
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Fig. 16 — Propagation curves for the first 16 modes for unity aspect ratio and 
An r -* 0. 



field lines reverse direction across the origin are no longer degenerate, 
but those whose field lines do not reverse still are degenerate. For all 
degeneracies to be split, there must exist a finite index difference as well 
as an aspect ratio other than unity. Figure 19 illustrates one such case. 
The effect of index difference on the degenerate principal modes for 
unity aspect ratio is examined in Fig. 20. The curve shows both a low 
and high index difference limit. In the range of interest for optical 




Fig. 17 — Propagation curves for the first 12 modes for a/b = 2 and An r — * 0. 
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Fig. 18 — Propagation curves for the first 16 modes for unity aspect ratio and 
An r = 0.5. 

circuits (0 — 0.1) the vanishing difference curve is an excellent ap- 
proximation. The greatest changes occur in the 0.1 — 10 range, which 
is the range of interest for some microwave problems. 

Figure 21 presents the computed results for the effect of index changes 
on the principal modes for an aspect ratio of 2. The effect is much 
stronger on the E\ x mode than the E x u mode. In fact, the effect on the 
E x u mode is comparatively small, except near cutoff. 

The effect of aspect ratio on the principal modes is demonstrated for 




Fig. 19 — Propagation curves for the first 12 modes for a/b = 2 and An, = 0.5. 
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Fig. 20 — E^ and #*, mode propagation curves for several values of An r with 
unity aspect ratio. 

vanishing index difference in Fig. 22. The curve for infinite aspect ratio 
was obtained from the exact analysis of the slab case. 14 

IV. CONCLUSIONS 

The results of the computations show that the circular harmonic 
method for analyzing rectangular dielectric waveguides gives excel- 




Fig. 21 — E\ x and E\ x mode propagation curves for several values of An r with 
o/6 = 2. 
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Fig. 22 — E" n and J5*, mode propagation curves for several values of a/b with 
An r -> 0. 

lent results for waveguides of moderate aspect ratio. The convergence 
of the computed results was rapid and the results are in agreement 
with those of Marcatili's in the regions where his approximations ap- 
ply. Furthermore, the results compare very well with Schlosser's 
curves for the principal mode. 

Comparison of the results presented here with Marcatili's show that 
the two methods give values of the normalized propagation constant, 
(P 2 , which are within a few percent for (P 2 > 0.5. Thus for (P 2 in this 
range his method is to be preferred since the calculations required are 
much simpler. However, for (P 2 < 0.5, and when it is desired to dif- 
ferentiate between modes for some of the near degenerate cases, an- 
other method must be used. 

The circular harmonic analysis is attractive for small (P 2 because of 
the nature of the matching boundary. For large refractive index dif- 
ference and moderate (P 2 both the method presented here and the one 
presented by Scholosser can be used. 
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